home *** CD-ROM | disk | FTP | other *** search
/ United Public Domain Gold 2 / United Public Domain Gold 2.iso / utilities / pu015.dms / pu015.adf / StarChart / Source / starcalc.c < prev    next >
C/C++ Source or Header  |  1989-09-24  |  14KB  |  510 lines
  1. /*=========================================================================
  2.   StarCalc.c -- This module contains the fixed radix arithmetic to speed up 
  3.                 the calculation of star positions. Floating point numbers
  4.             are converted to the fixed radix format by multiplying by
  5.             10000, and the result is stored in a LONG integer. The
  6.             routine sacrifices space (and some accuracy) in the interests
  7.             of speed (the old tradeoff). The FastSin, FastCos, and
  8.             FastTan functions rely on a table of 900 precalculated
  9.         sines and tangents (in TrigTab.h) that are in the fixed
  10.         format. 
  11.   
  12.   Credits for Star Chart:
  13.        Robert L. Hill of the Orange County, CA. Amiga Friends User Group
  14.                       wrote the original version of StarChart in AmigaBasic
  15.                       The star data and many of the main functions of this
  16.                       version are derived from that program.
  17.  
  18.        Ray R. Larson  wrote the c version 1.0 of StarChart, 'intuitionizing'
  19.                       and enhancing the speed and functions of the original.
  20.  
  21.   Copyright (c) 1986 by Ray R. Larson
  22.   
  23.   This program may be freely distributed and copied, but may not be sold
  24.   without the permission of the author. If you modify or enhance it, 
  25.   please include the above credits (and please send me a copy!).
  26.  
  27. Ray R. Larson
  28. 6425 Central Ave. #304
  29. El Cerrito, CA 94530
  30.  
  31. BitNet  LARSON@UCBCMSA
  32. Well    rrl
  33. CServe  70446,766
  34. =========================================================================*/
  35. /*------------Header file for all of the standard stuff----*/ 
  36. /*-------------plus definitions of global structures-------*/
  37. #include "star.h" 
  38.  
  39. /*------------Header file for sin/cos/tan/cotan table -----*/
  40. #include "TrigTab.h"
  41.  
  42. SHORT FirstStar, LastStar;
  43.  
  44. extern struct ParamStruct Parms;
  45. extern struct Coord coords[]; /* x and y screen coordinates for stars */
  46. extern struct star_rec StarTable[];
  47. extern FLOAT P1,P12,P15; /* pi based values */
  48.  
  49.  
  50. /* this routine fetches the sine values for angles specified in tenths of */
  51. /* a degree (i.e. 451 is 45.1 degrees)                                    */
  52. LONG FastSin(tenthdegs)
  53. LONG tenthdegs;
  54. {
  55.  
  56.    register LONG deg, sign;
  57.  
  58.    if (tenthdegs == 0L) return(0L);
  59.  
  60.    if (tenthdegs < 0L) 
  61.       { sign = -1L;
  62.         deg = -tenthdegs;
  63.       }
  64.    else 
  65.       { sign = 1L;
  66.         deg = tenthdegs;
  67.       }
  68.    
  69.       
  70.    deg = deg % 3600L;
  71.    if (deg > 1800L)
  72.       { sign =  (sign == -1L)? 1L : -1L;
  73.         deg = deg - 1800L;
  74.       }
  75.     if (deg > 900L) deg = (900L-deg)+900L;
  76.     
  77.     return(SinTanTab[deg].sine * sign);
  78. }
  79.  
  80. /* FastCos just shifts the sine curve and calls FastSin */
  81. LONG FastCos(tenthdegs)
  82. LONG tenthdegs;
  83. {
  84.    return(FastSin(tenthdegs + 900L));
  85. }
  86.  
  87. /*  FastTan looks up the tangent  */
  88. LONG FastTan(tenthdegs)
  89. LONG tenthdegs;
  90. {
  91.    register LONG deg, sign;
  92.  
  93.    if (tenthdegs == 0L) return(0L);
  94.  
  95.    deg = (tenthdegs < 0L) ? -tenthdegs : tenthdegs;
  96.    deg = deg % 1800L;
  97.  
  98.    if (deg < 900L) return(SinTanTab[deg].tangent);
  99.    else 
  100.        return(SinTanTab[(900L - deg) + 900L].tangent * -1L);
  101. }
  102.  
  103. /* convert float to fixed format */
  104. LONG ftofix(x)
  105. FLOAT x;
  106. { return((LONG)(x * 10000.0));
  107. }
  108.  
  109. /* multiple two fixed radix numbers */
  110. LONG fixmult(x,y)
  111. LONG x, y;
  112. {
  113.     long temp1, temp2;
  114.     temp1 = (x * (y % 10000))/10000;
  115.     temp2 = x * (y/10000);
  116.     
  117.     return (temp1+temp2);
  118. }
  119.  
  120. /* divide two fixed radix numbers , lose a little precision here...*/
  121. LONG fixdiv(x,y)
  122. LONG x, y;
  123. {
  124.     long temp;
  125.     if (x == 0L) return (0L);
  126.     if (x == y) return (10000); /* i.e. 1 */
  127.     /* shift x over by 2 digits */
  128.     
  129.     temp = x * 100;
  130.     return (temp/y * 100);
  131. }
  132.  
  133. /***************************************************************************
  134.  * FastCalc - Calculate the positions of the visible stars given the
  135.  *  current parms and plot them on the screen. Save positions in coords.
  136.  *  This version uses fixed radix arithmetic
  137.  ***************************************************************************/
  138. FCalc1()
  139. {
  140.   LONG xDeg,yDeg,Xpos,Ypos, LatDeg, LatPos, LatCOS, LatSIN, F, LST;
  141.   SHORT i;
  142.   LONG yeardif, RA, DC;
  143.   BOOL visible;
  144.   struct star_rec *ST;
  145.   struct Coord *c;
  146.  
  147.   /* set some (variable) constants */
  148.   yeardif = Parms.Year - 1950;
  149.   LST = ftofix(Parms.Lst);
  150.   LatDeg = (LONG)(Parms.Latitude * 10.0);
  151.   LatPos = 1600000 - fixmult(17800, ftofix(Parms.Latitude));
  152.   LatCOS = FastCos(LatDeg);
  153.   LatSIN = -FastSin(LatDeg);
  154.   F = 891268;  /* ie 140/(PI/2) * 10000 */
  155.  
  156.   FirstStar = 0;
  157.   
  158.   /* for a little extra speed we will use pointers rather than array indexes */
  159.   ST = &StarTable[1];
  160.   c = &coords[1];
  161.  
  162.   for (i = 1; i <= NumStars; i++)
  163.     { 
  164.      if (yeardif) {
  165.        RA = ftofix(ST->Ra);
  166.        DC = ftofix(ST->Dc);
  167.        /* convert right Asc. and declination of star to tenth Degrees */
  168.        xDeg = (RA * 15)/1000;
  169.        yDeg = DC/1000;
  170.        
  171.      
  172.        /* correct for year - uses mathffp routines*/
  173.        Xpos = RA + 
  174.              (30700 + fixmult(fixmult(13400, FastSin(xDeg)), tan(xDeg)))
  175.           * yeardif / 3600;
  176.  
  177.       
  178.        Ypos = DC + fixmult(200000, FastCos(yDeg)) * yeardif / 3600;
  179.       }
  180.    else { /* year is 1950 - basis for the star table */
  181.          Xpos = RA;
  182.          Ypos = DC;
  183.         }
  184.   
  185.    /* modify the Xpos and Ypos for the date, time, location and horizon */   
  186.    if(Parms.Horizon == 0) 
  187.      visible = FastNorth(&Xpos,&Ypos,LatDeg,LatPos,LatCOS,LatSIN,F,LST);
  188.    else 
  189.      visible = FastSouth(&Xpos,&Ypos,LatDeg,LatCOS,LatSIN,LST);
  190.   
  191.    /* if the star is visible from the current location - save coords and   */
  192.    /* plot it on the display (not any more - actual plot is done in redisp */
  193.    if(visible)
  194.      { /* set the display coordinates */
  195.        c->x = CHARTLEFT + (((2 * Xpos) + 5000)/10000);
  196.        c->y = (CHARTTOP+1L) + ((Ypos + 5000)/10000);
  197.  
  198.        /* set the first and last coords filled in */
  199.        if (FirstStar) LastStar = i;
  200.        else FirstStar = i;
  201.      }
  202.    else 
  203.      { c->x = 0L;
  204.        c->y = 0L;
  205.      }
  206.      
  207.    /* increment the pointers */
  208.    ST++;
  209.    c++;  
  210.    } /* end of for loop */
  211. } /* end of DisplayStars */
  212.  
  213.  
  214. /***************************************************************************
  215.  * FastNorth - calculate star position for Northern Sky
  216.  *  returns 1 if star is visible, zero if not visible from current location
  217.  ***************************************************************************/
  218. BOOL FastNorth (xposit,yposit,lat,LatPos,latcos,latsin,F,LST)
  219. LONG *xposit, *yposit,lat,LatPos,latcos,latsin,F,LST;
  220. {
  221.   LONG DecDeg,HourDeg, Temp1, TempX, TempY;
  222.   register LONG xpos, ypos;
  223.  
  224.   xpos = *xposit;
  225.   ypos = *yposit;
  226.  
  227.   if (ypos < 0) return(0);
  228.   
  229.   /* convert declination, latitude, and local siderial time to Degrees */
  230.   DecDeg = ypos/1000;
  231.   HourDeg = (LST - xpos)*15 / 1000;
  232.   /* eliminate stars below the horizon line */
  233.   if(( fixmult(latcos,fixmult(FastCos(DecDeg), FastCos(HourDeg))))
  234.      <( fixmult(latsin,FastSin(DecDeg))))
  235.       return(0);
  236.  
  237.   xpos = xpos - LST;
  238.   
  239.   /* Map to scale for the display */
  240.   xpos = (xpos * 15) - 15708;
  241.  
  242.   Temp1 =  fixmult(F, (15708 -  abs(ypos)));
  243.   TempX = fixmult(Temp1,FastCos((xpos/1000))) + 1400000;
  244.   
  245.   TempY = fixmult(Temp1, FastSin((xpos/1000))) + LatPos;
  246.   
  247.   xpos = TempX;
  248.   ypos = TempY;
  249.  
  250.   if((xpos > 2790000) || (xpos < 0)) return(0);
  251.   if((ypos > 1590000) || (ypos < 0)) return(0);
  252.   
  253.   /* if we got to here it must be visible */
  254.   *xposit = xpos;
  255.   *yposit = ypos;
  256.   return(1);
  257. } /* end of FastNorth */
  258.  
  259. /***************************************************************************
  260.  * FastSouth - calculate star position for Southern Sky
  261.  *  returns 1 if star is visible, zero if not visible from current location
  262.  ***************************************************************************/
  263. BOOL FastSouth (xposit,yposit,latdeg,latcos,latsin,LST)
  264. LONG *xposit, *yposit,latdeg,latcos,latsin,LST;
  265. {
  266.   LONG DecDeg,LatDeg,HourDeg, LatPos, Temp1;
  267.   register LONG xpos, ypos;
  268.  
  269.   xpos = *xposit;
  270.   ypos = *yposit;
  271.  
  272.   if (ypos > (latdeg*1000)) return(0);
  273.  
  274.   if (ypos <  ((latdeg*1000) - 900000)) return(0);
  275.  
  276.   
  277.   /* convert declination, latitude, and local siderial time to Degrees */
  278.   DecDeg = ypos/1000;
  279.   HourDeg = (LST - xpos)*15 / 1000;
  280.   /* eliminate stars below the horizon line */
  281.   if(( fixmult(latcos,fixmult(FastCos(DecDeg), FastCos(HourDeg))))
  282.      <( fixmult(latsin,FastSin(DecDeg)))) return(0);
  283.  
  284.   xpos =  LST - xpos;
  285.   if(xpos < -120000) xpos += 240000;
  286.   if(xpos > 120000) xpos -= 240000;
  287.   xpos = 1400000 + fixmult(233300 ,xpos);
  288.   
  289.   if((xpos > 2790000) || (xpos < 0)) return(0);
  290.  
  291.   ypos = fixmult(17800, (latdeg*1000) - ypos);
  292.   
  293.   if((ypos > 1590000) || (ypos < 0)) return(0);
  294.   
  295.   /* if we got to here it must be visible */
  296.   *xposit = xpos;
  297.   *yposit = ypos;
  298.   return(1);
  299. } /* end of FastSouth */
  300.  
  301.  
  302. FLOAT FSin(x)
  303. FLOAT x;
  304. {
  305.    LONG tenthdegs;
  306.    
  307.    
  308.    tenthdegs = (LONG) (x * 57.2958)* 10.0; /* convert to degrees */
  309.    return(((FLOAT)FastSin(tenthdegs)/ 10000.0));
  310. }
  311.  
  312.  
  313. FLOAT FCos(x)
  314. FLOAT x;
  315. {
  316.    LONG tenthdegs;
  317.    
  318.    
  319.    tenthdegs = (LONG) (x * 57.2958)* 10.0; /* convert to degrees */
  320.    return(((FLOAT)FastCos(tenthdegs)/ 10000.0));
  321. }
  322.  
  323. FLOAT FTan(x)
  324. FLOAT x;
  325. {
  326.    LONG tenthdegs;
  327.    
  328.    
  329.    tenthdegs = (LONG) (x * 57.2958)* 10.0; /* convert to degrees */
  330.    return(((FLOAT)FastTan(tenthdegs)/ 10000.0));
  331. }
  332.  
  333.   
  334.  
  335. /***************************************************************************
  336.  * CalcStars - Calculate the positions of the visible stars given the
  337.  *  current parms and plot them on the screen. Save positions in coords.
  338.  ***************************************************************************/
  339. FastCalc()
  340. {
  341.   FLOAT xRad,yRad,Xpos,Ypos, LatRad, LatPos, LatCOS, LatSIN, F;
  342.   SHORT i;
  343.   FLOAT yeardif, RA, DC;
  344.   BOOL visible;
  345.   struct star_rec *ST;
  346.   struct Coord *c;
  347.  
  348.   /* set some (variable) constants */
  349.   yeardif = (FLOAT)(Parms.Year - 1950);
  350.   LatRad = Parms.Latitude * P1;
  351.   LatPos = 160.0 - (1.78 * Parms.Latitude);
  352.   LatCOS = FCos(LatRad);
  353.   LatSIN =  -FSin(LatRad);
  354.   F = 140.0 / 1.5708;
  355.  
  356.   /* for a little extra speed we will use pointers rather than array indexes */
  357.   ST = &StarTable[1];
  358.   c = &coords[1];
  359.  
  360.   for (i = 1; i <= NumStars; i++)
  361.     { 
  362.      if (((Parms.Horizon == 0)&&(ST->Dc < 0.0)) || /* north view */
  363.          ((Parms.Horizon == 1)&&
  364.           ((ST->Dc > (FLOAT)Parms.Latitude)|| 
  365.            (ST->Dc < (FLOAT)(Parms.Latitude - 90))) /* south view */
  366.           ))
  367.         { /* don't bother calculating if it will never be visible given the */
  368.           /* current position and orientation.                    */
  369.           c->x = 0L;
  370.           c->y = 0L;
  371.           ST++;
  372.           c++;  
  373.           continue;
  374.          }
  375.          
  376.      if (yeardif) {
  377.        RA = ST->Ra;
  378.        DC = ST->Dc;
  379.        /* convert right Asc. and declination of star to radians */
  380.        xRad = RA * P12;
  381.        yRad = DC * P1;
  382.        
  383.      
  384.        /* correct for year - uses mathffp routines*/
  385.        Xpos = RA + 
  386.              (3.07 + 1.34 * FSin(xRad) * FTan(xRad))
  387.           * yeardif / 3600.0;
  388.  
  389.       
  390.        Ypos = DC + 20.0 * FCos(yRad) * yeardif / 3600.0;
  391.       }
  392.    else { /* year is 1950 - basis for the star table */
  393.          Xpos = RA;
  394.          Ypos = RA;
  395.         }
  396.   
  397.    /* modify the Xpos and Ypos for the date, time, location and horizon */   
  398.    if(Parms.Horizon == 0) 
  399.      visible = NorthVF(&Xpos,&Ypos,LatRad,LatPos,LatCOS,LatSIN,F);
  400.    else 
  401.      visible = SouthVF(&Xpos,&Ypos,LatCOS,LatSIN);
  402.   
  403.    /* if the star is visible from the current location - save coords and   */
  404.    /* plot it on the display (not any more - actual plot is done in redisp */
  405.    if(visible)
  406.      { /* set the display coordinates */
  407.        c->x = CHARTLEFT + (LONG)((2.0 * Xpos) + 0.5);
  408.        c->y = (CHARTTOP+1L) + (LONG)(Ypos + 0.5);
  409.      }
  410.    else 
  411.      { c->x = 0L;
  412.        c->y = 0L;
  413.      }
  414.      
  415.    /* increment the pointers */
  416.    ST++;
  417.    c++;  
  418.    } /* end of for loop */
  419. } /* end of DisplayStars */
  420.  
  421.  
  422. /***************************************************************************
  423.  * NorthView - calculate star position for Northern Sky
  424.  *  returns 1 if star is visible, zero if not visible from current location
  425.  ***************************************************************************/
  426. BOOL NorthVF(xposit,yposit,lat,LatPos,latcos,latsin,F)
  427. FLOAT *xposit, *yposit,lat,LatPos,latcos,latsin,F;
  428. {
  429.   FLOAT DecRad,HourRad, Temp1, TempX, TempY;
  430.   register FLOAT xpos, ypos;
  431.  
  432.   xpos = *xposit;
  433.   ypos = *yposit;
  434.  
  435.   if (ypos < 0.0) return(0);
  436.   
  437.   /* convert declination, latitude, and local siderial time to radians */
  438.   DecRad = ypos * P1;
  439.   HourRad = P12 * (Parms.Lst - xpos);
  440.   
  441.   /* eliminate stars below the horizon line */
  442.   if(( latcos * FCos(DecRad) * FCos(HourRad))
  443.      <( latsin * FSin(DecRad))) return(0);
  444.  
  445.   xpos = xpos - Parms.Lst;
  446.   ypos = ypos * P1;
  447.   
  448.   /* Map to scale for the display */
  449.   xpos = (xpos * P15) - 1.5708;
  450.  
  451.   Temp1 =  F * (1.5708 -  abs(ypos));
  452.   TempX = Temp1 * FCos(xpos) + 140.0;
  453.   
  454.   TempY = Temp1 * FSin(xpos) + LatPos;
  455.   
  456.   xpos = TempX;
  457.   ypos = TempY;
  458.  
  459.   if((xpos > 279.0) || (xpos < 0.0)) return(0);
  460.   if((ypos > 159.0) || (ypos < 0.0)) return(0);
  461.   
  462.   /* if we got to here it must be visible */
  463.   *xposit = xpos;
  464.   *yposit = ypos;
  465.   return(1);
  466. } /* end of NorthView */
  467.  
  468. /***************************************************************************
  469.  * SouthView - calculate star position for Southern Sky
  470.  *  returns 1 if star is visible, zero if not visible from current location
  471.  ***************************************************************************/
  472. BOOL SouthVF (xposit,yposit,latcos,latsin)
  473. FLOAT *xposit, *yposit,latcos,latsin;
  474. {
  475.   FLOAT DecRad,LatRad,HourRad, LatPos, Temp1;
  476.   register FLOAT xpos, ypos;
  477.  
  478.   xpos = *xposit;
  479.   ypos = *yposit;
  480.  
  481.   if (ypos > (FLOAT)Parms.Latitude) return(0);
  482.  
  483.   if (ypos < (FLOAT)(Parms.Latitude - 90)) return(0);
  484.  
  485.   
  486.   /* convert declination, latitude, and local siderial time to radians */
  487.   DecRad = ypos * P1;
  488.   HourRad = P12 * (Parms.Lst - xpos);
  489.   
  490.   /* eliminate stars below the horizon line */
  491.   if(( latcos * FCos(DecRad) * FCos(HourRad))
  492.      <( latsin * FSin(DecRad))) return(0);
  493.  
  494.   xpos =  Parms.Lst - xpos;
  495.   if(xpos < -12.0) xpos += 24.0;
  496.   if(xpos > 12.0) xpos -= 24.0;
  497.   xpos = 140.0 + (23.33 * xpos);
  498.   
  499.   if((xpos > 279) || (xpos < 0)) return(0);
  500.  
  501.   ypos = 1.78 * (Parms.Latitude - ypos);
  502.   
  503.   if((ypos > 159.0) || (ypos < 0.0)) return(0);
  504.   
  505.   /* if we got to here it must be visible */
  506.   *xposit = xpos;
  507.   *yposit = ypos;
  508.   return(1);
  509. } /* end of SouthView */
  510.